< 



High-Precision Numerical Determination of 
Eigenvalues for a Double- Well Potential Related to 
in" the Zinn- Justin Conjecture 

o 

(N . H. A. Alhendi 1 and E. I. Lashin 1 ' 2 

1 Department of physics and Astronomy, College of Science, 
King Saud University, Riyadh, Saudi Arabia 
2 Department of Physics, Faculty of Science, 

Ain Shams University, Cairo, Egypt 
Email: alhendi@ksu.edu. sa, lashin@ksu.edu. sa 

; February 1, 2008 



(N 

A numerical method of high precision is used to calculate the energy eigenvalues and 
eigenfunctions for a symmetric double-well potential. The method is based on enclosing the 
system within two infinite walls with a large but finite separation and developing a power 
series solution for the Schrodinger equation. The obtained numerical results are compared 
with those obtained on the basis of the Zinn-Justin conjecture and found to be in an excellent 
c"j ■ agreement. 
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k> ■ 1 Introduction 

Quantum mechanical tunneling through finite barriers is a well established phenomenon in theory 
and application. The symmetric double well potential is one of the many examples exhibiting 
this phenomenon. In this case, the energy splitting generated by tunneling can be estimated 
with the help of the well-known semi-classical WKB approximation and instanton techniques 
(see for example [1]). However, to calculate this splitting accurately, one needs an effective 
method of high precision. 

In a series of papers, Zinn-Justin [2] developed a conjecture (to be termed 'the Zinn-Justin 
conjecture') to determine the energy levels of a quantum Hamiltonian H, in cases where the 
potential has degenerate minima. This conjecture takes the form of the generalized Bohr- 
Sommerfeld quantization formulae. It has been applied, among other potentials, to the case of 
the symmetric double well. In this case the Hamiltonian is 

H = ~^-?2+- V (l), where V {q) = \q 2 {I - q) 2 . (1) 
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It is obvious that this Hamltonian is invariant under the transformation (q — > 1 — q). The energy 
eigenvalues for this potential can be obtained by finding a solution to the Zinn-Justin conjecture 
equation: 

1/1 \ / 2\ D ( E >^ 

,T --D(E,g))—) exp[-A(E,g)/2]=±i. (2) 
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The function D(E,g) has a perturbative expansion in powers of g, of which the first few terms 



are 



D(E, g)=E + g (sE 2 + \)+ 9 2 (<S5E 3 + ^ + 0(g 3 ). (3) 

The other function A(E, g) receives contributions from the instanton expansion in the path 
integral and its first few terms are 

A(E, 9) = y g +9 (l7E 2 + + g 2 ^27E 3 + + 0(g 3 ). (4) 

The energy E^,± can be extracted from equation (2) by expanding in powers of g and in the 
two quantities 

A( 9 )=.n(-H) and «,) = f^I^M. (5) 
The complete semi-classical expansion of Epj ± has the form [4] 

oo oo /oxTVn/ „-l/6p\ n n-l oo ,„ r , , 

^w = E* + E EM-Vs)) fc E e 2 V. (6) 

Z=0 n=l Vi// V / fc= ;=o 

The coefficients e relevant to the numerical calculation have been explicitly calculated in [3]. 
The number N is the unperturbed quantum number which corresponds to 

E±, N (g) = N +1/2 + 0(g). (7) 

A detailed exposition of the above equations can be found in [4] . 

In [3] , numerical calculations have been carried out and led to the energy eigenvalues for the 
ground and the first excited states respectively, for g = 0.001, 

£ ,+(0.001) = 0.49899 54548 62109 17168 91308 39481 92163 68209 47240 20809 

66532 93278 69722 01391 15135 28505 38294 45798 45759 95999 

06739 55175 84722 67802 81306 96906 01325 25943 77289 94365 

88255 24440 17437 12789 27978 99793 , (8) 

Eq- (0.001) = 0.49899 54548 62109 17168 91308 39481 92163 68209 47240 20809 

66532 93278 69722 01391 29839 92959 55803 70812 27749 92448 

48259 36743 64757 68328 84835 35511 34663 06309 82331 51885 

23308 08622 84780 52722 10103 67282 . (9) 

The above numerical results have been obtained by lattice extrapolation using a modified 
Richardson algorithm [3]. 

This tiny difference encourages us to seek for an independent but simple and direct method, 
which allows us to obtain the energy eigenvalues for the potential in equation (1) and compare 
them with the above numerical results. In addition, the present method allows us to obtain 
an accurate description for the corresponding wavefunctions. This method has been previously 
applied to various potential functions with and without degenerate minima, leading to results 
with high accuracy [5]. 

The method, as will be described in the next two sections, is based on power series solution 
of the Schrodinger equation in a finite range. It has appeared from time to time in the literature 
[6, 7, 8], but has not been developed to its maximum efficiency. We shall show that, by using 
the computer algebra systems (for example Mathemtica) which can deal with exact numbers, 
the accuracy of the method can be substantially improved. 

In the following section, for illustrative purpose, we explain our method by applying it to the 
well-known exactly solvable harmonic oscillator potential and then extend it to the symmetric 
double well. 



2 Calculations and Results 



In this section we, first, consider the well-known exactly solvable harmonic oscillator. In this 
case, the Schrodinger equation reads (h = l,m = 1) 



2dq 



2 +E-V(q) 



= 0, (10) 



where 



V{q) = \ Q 2 . (11) 
The exact energy eigenvalues and the corresponding eigenfunctions are 

En = ( N +l)> JV = 0,1,2,---, 

V N (q) = 2-T (iV!)-5 exp (-^J i^), 

(12) 

where HN^q) are the Hermite polynomials. 

For the harmonic oscillator confined between two infinite walls at q = ± L, we develop a 
power series solution in the form 

oo 

%) = E a *?' ( 13 ) 

i=0 

Substituting in equation (10), one gets the following recursion relation: 

di = a% ~' 4 , - ^ a% ~ 2 i + o, 1 and aj = when i < 0. (14) 
i(t-l) 

The symmetry of the potential implies that we have two types of solutions: the even solutions 
obtained by imposing (ignoring normalization) oo = l,ai = and the odd ones by imposing 
ao = 0, a\ = 1. The energy eigenvalues are then obtained from the condition ty(E,L) = for 
both cases. 

For numerical calculations, we approximate the power series in equation (13) with a truncated 
one having a finite number of terms *$>i(E,q), where / is the number of non- vanishing terms. 
The boundary condition for a specific value of L corresponds to ^i(E, L) = 0. To get the zeros 
of ^j(E,L) with respect to E, we first plot a graph for tyj(L,E) as a function of E to locate 
where ^i(L,E) changes sign. We then can use two nearby points containing one single root 
as the initial iteration for the 'bisection method' to find the zeros. In doing this we have used 
Mathematica package version 3 and also have relied extensively on its ability to manipulate 
exact numbers. The stability of the numerical results to a certain degree of accuracy is checked, 
for a particular L, by increasing / till the obtained value of E stays fixed. 

In table 1 we present the calculated energies for the ground and the first three excited states 
for the bounded harmonic oscillator as compared to the exact results of the unbounded one. 



L N E 



N 



^exact ~ (jV ^ry 



250 8 0.50000000000000000000000000 

1 1.500000000000000000000000 

2 2.5000000000000000000000 

3 3.500000000000000000000 



Table 1: The calculated first four energy levels for the bounded harmonic oscillator compared to the 
unbounded one. (2L) is the width of the well and I refers to the number of the non-vanishing terms in 
the truncated series of the wavefunction. 

The present method, as can be seen from table 1, reproduces for a large value of L, the 
exact ones even for a moderate number of non-vanishing terms in the truncated series of the 
wavefunction. Moreover, one can get an accurate description for the wavefunctions shown in 
figure 1 which can not be distinguished from the exact ones when drawn within the same interval 
\q\ < L = 8. 
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Figure 1: The normalized ground (left) and first excited (right) state wavefunctions for the bounded 
harmonic oscillator for L — 8 . 

Now we apply the above-explained method to the double well potential in equation (1). For 
our convenience, we use the substitution q — > q + ^, so the potential in equation (1) now takes 
the form 

V(q) = \ (q + 1/2) 2 (q - 1/2) 2 . (15) 

This form of the potential has now inversion symmetry (q — > —q) which is suitable for our 
calculation. It should be evident that rewriting the potential in this form doesn't affect the 
eigenvalues of the Hamiltonian in equation (1). As explained above, for this potential we again 
use the power series expansion of the wavefunction in the finite range. The Schrodinger equation, 
for the potential V(q) in equation (15), is 

2 aq z g 
In substituting the power 

i 



V(q) =0, —L < q < L. (16) 

series expansion, 

(17) 



in equation (16), one gets the following recurrence formula for the expansion coefficients, af. 

-6 — + Je a i-2 — Edi-2 



J_ 

2ff 



i(i - 1) 

For L = 3, the obtained eigenvalues are 
£ ,+ (0.001) 



i 0,1 and Oj = when i < 0. 



£b,_ (0.001) 



0.49899 54548 62109 17168 91308 39481 92163 68209 47240 20809 
66532 93278 69722 01391 15135 28505 38294 45798 45759 95999 
06739 55175 84722 67802 81306 96906 01325 25943 77289 94365 
88255 24440 17437 12789 27978 99793 98922 00536 06978 04138 
65255 73028 37723 50241 67171 . 

0.49899 54548 62109 17168 91308 39481 92163 68209 47240 20809 
66532 93278 69722 01391 29839 92959 55803 70812 27749 92448 
48259 36743 64757 68328 84835 35511 34663 06309 82331 51885 
23308 08622 84780 52722 10103 67282 72047 61340 01672 24803 
65523 52410 13798 16304 58360. 



(18) 



(19) 



(20) 



These values agree with the ones obtained from the numerical calculations based on the Zinn- 
Justin conjecture. In figure 2 we present the ground and the first excited state wavefunctions 
for the bounded double-well potential for g = 1/1000, I = 4600 and L = 1. 
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Figure 2: The normalized ground (left) and first excited (right) state wavefunctions for the bounded 
double-well potential with g = 1/1000 . 



3 Discussion 

It is important to note the following generic remarks. First, a reason for the capability of the 
present method is that for a bound state, the wavefunction is spatially localized which means 
that the probability density (l^) 2 ) has appreciable values in a finite region of space behind which 
the probability density tends rapidly to zero. Thus, to a good approximation, it is, therefore, 
reasonable to consider the corresponding problem in a finite interval, with a suitable width, 
bounded by two infinite walls. The criteria for a suitable value of L can be quantitatively given 
by the condition E « V(L). Second, from the WKB approximation, it can be made plausible 
that the zeros of ^(E, L) provide upper bounds for the energy eigenvalues while the zeros of the 
derivative-with respect to q— *f>'(E', L) provide the lower ones; the same finding can be proved 
in a rigorous way as shown in [9]. Thus, by matching the digits of the two zeros, one can get an 



accurate energy eigenvalue up to the number of coincident digits. As an example for the ground 
state of the bounded harmonic oscillator, with L = 8 and / = 250, one gets 

E = 0.5000000000000000000000000014362707054755765903756598 

26757972824824621785332078167891514939744867648, 
E' = 0.4999999999999999999999999985405543573278682092744652 

58622103903146216005437303539479001558808137418. 

The corresponding wavefunctions and their slopes are 

V(E, 8) = 4.8 x 10~ 49 , 8) = 8.1 x 10~ 48 , 

V(E', 8) = 2.6 x 10~ 14 , 8) = -2 x 10~ 13 . 

(22) 

After matching the digits of the two numbers in equation (21), one gets the ground-state energy 
accurate up to 25 digits as shown in table 1. The remaining eigenvalues are obtained by the same 
procedure. However, one should pay attention that this accuracy is expected to be less than the 
accuracy of the bisection method. In this method, the accuracy estimation is e = (c — a)/2 n 
where n, here, is the number of iteration, and c and a are two points enclosing only one root. 
In our case, we have taken for the ground state n = 200, c = yjj and a = 4/10 , giving 
e = 1.2 x 10 -61 . Finally according to the WKB approximation, the wavefunction behaves for 
large q in the inaccessible region as 

^WKB(q) oc — — ~y exp (- f J(V(q')-E) dq') , (23) 

(V(q) — E)i \ Jqt J 

where qt is a turning point just left to the inaccessible region. The value of ^wkb{q = 8) is 6.5 x 
10~ 14 , while for the truncated series solution of equation (13) it has the value 4.8 x 10 -49 as 
given in equation (22). The reason for this huge difference is that the series solution is valid and 
convergent as long as q is finite [10]. In addition to this, the energy eigenvalues as extracted from 
the zeros of ty(E,L) (for suitable L) result in a delicate cancellation between terms of opposite 
signs in the power series solution. 

One may suspect that using a series solution in the form 

i&(x) = exp (— bx 2 ) V] cij x J , 

3 

may help improving the rate of convergence for the obtained eigenvalues. In contrast, one needs 
more terms in the series expansion to achieve the same level of accuracy obtained by the series 
solution of the form given in equation (13). The reason behind this stems from the fact that 
any finite truncation for the series in the form given in equation (24) always decays, due to the 
exponential factor, as q becomes large making the determination of the energy eigenvalues less 
reliable, especially when the parameter b is large. As an example, when b = \ we can achieve the 
same accuracy reported in table 1 with the same number of non-vanishing terms in the truncated 
series, while for b = 8 we need 600 non-vanishing terms to achieve the same accuracy. Thus, the 
best thing which can be done is to work with the parameter b having zero value. However, it 
should be kept in mind that both series in equation (13) and equation (24) are equivalent but 
only in the infinite sum limit. 

We also study the effect of the parameter b in the case of the double-well potential given by 

V(x) = -10 x 2 + x 4 , (in units ft=l,m=-J. (25) 



(21) 



(24) 



As an example, when we work with the precision 100 digits, then we find for 6 = 0, I = 
750 and L = 8, that the ground state energy has the value (accurate up to 69 digits) 

E Q = -20.63357 67029 47799 14995 85548 37431 

50876 53159 46057 73551 39057 10311 42892 92. (26) 

To achieve the same accurate energy determination for 6 = 10, we find that it is possible to 
use 500 terms which is not considerably less than the case of 6 = 0. However, this comes with 
the high cost of working with precision 300 digits. Working with such a high precision renders 
the calculation slow. At intermediate values of b like 2, 3 and 4, we can use less terms but with 
high precision as shown in table 2. According to our numerical investigations for the case of the 
double well, in the finite range, the choice b = is the best compromise between the number of 
terms used and the degree of precision to get a more efficient calculation. 
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750 


750 


750 


500 


500 


500 


500 


500 


Precision 


100 


100 


150 


150 


200 


200 


200 


300 



Table 2: Precision versus I and the parameter b 

It is important to point out that in dealing with low accuracy results (like nine digits), one 
cannot decide which is better, to work with or without the parameter b. Furthermore, employing 
the method in a non-efficient way may lead to wrong conclusions as in [8] , where it is emphasized 
that setting a non-vanishing value for the parameter b greatly reduces the number of terms 
used. To clarify these points, we obtain for the potential given by equation (25) the four first 
energy levels (E = -20.6335767, E x = -20.6334568, E 3 = -12.3795437, E 4 = -12.3756738) 
accurate up to 10 digits as presented in [8]; our results (using L = 4.2) are summarized in table 3. 
It is evident from table 3 that one can not say it is a big advantage to use 90 terms (for 6=2) 
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125 


100 


90 


90 


90 


90 


90 


200 



Table 3: The parameter b versus I (number of non- vanishing terms) 

rather than 125 terms (for 6 = 0). However, numerical studies clearly indicate that the situation 
becomes worse when 6 increases (for 6 = 10 we need 200 terms). Another clear example is the 
pure quartic potential (V(x) = x 4 ) for which we get, for 6 = 0, L = 3.5, and I = 75, low- 
energy eigenvalues (the first five) determined accurately up to nine digits while obtaining the 
same results for the choice 6 = 3 and I = 50. Furthermore, the tenth eigenvalue is determined 
accurately up to 9 digits, for L = 3.9, using / = 75 for 6 = 3, while / = 125 for 6 = 0. These 
findings are in contradiction with what has been claimed in [8], where it was mentioned that one 
should use about 2000 terms in the power series to determine the energy for the choice 6 = 0. 
Similar findings occur for the potential V(x) = x 2 + a: 8 . In such a situation for L = 2.5, we can 
use 125 terms in the power series solution for 6 = and 75 terms for 6 = 5, while getting the 
same accurate results up to nine digits. 

The problem in the calculations found in [6, 8] comes from evaluating every term in the power 
series to a certain precision, and then summing the series which leads to an error accumulation, 



resulting in low-accuracy results despite using a large number of terms. In our approach, we sum 
all terms in the power series exactly, and then only in determining the roots (energy) from the 
condition ^i(E, L) = 0, do we resort to numerical calculation with a certain precision. Although 
the ability of the computer algebra system to deal with exact numbers was available from the 
early 1980s, it has not been used since then in such calculations. 

4 Conclusion 

In this paper we have presented an independent simple method leading to eigenvalues which 
agree well with the recently obtained numerical results based on the Zinn-Justin conjecture for 
the symmetric double-well potential. We have also included results with more significant digits 
than reported. It has been applied to some other potentials to illustrate its capability, and 
its precision has been compared with other calculations based on introducing an exponentially 
decaying factor (e~ bx ). Several subtle points related to its precision have also been discussed 
and clarified. The method we opted also enables us to get an accurate numerical determination 
of the corresponding wavefunctions. 
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